#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <ctype.h>
#include <libxml/xmlreader.h>
#include <libxml/xmlwriter.h>
#include <complex.h>
#include <fftw3.h>
#include <argp.h>
#include "RPFITS.h"
#include "cpgplot.h"

#define BUFSIZE 256

#define YES 1
#define NO  0

/* XML encoding type */
#define MY_ENCODING "ISO-8859-1"

/* RPFITS commands numerical definitions to make the code more readable */
#define JSTAT_OPENFILE -3
#define JSTAT_OPENFILE_READHEADER -2
#define JSTAT_READNEXTHEADER -1
#define JSTAT_READDATA 0
#define JSTAT_CLOSEFILE 1
#define JSTAT_SKIPTOEND 2

/* RPFITS return codes definitions to make the code more readable */
#define JSTAT_UNSUCCESSFUL -1
#define JSTAT_SUCCESSFUL 0
#define JSTAT_HEADERNOTDATA 1
#define JSTAT_ENDOFSCAN 2
#define JSTAT_ENDOFFILE 3
#define JSTAT_FGTABLE 4
#define JSTAT_ILLEGALDATA 5

/* the different colour tables we can use in PGPLOT */
#define PGPLOT_CT_GREYSCALE 1
#define PGPLOT_CT_RAINBOW 2
#define PGPLOT_CT_HEAT 3
#define PGPLOT_CT_WEIRDIRAF 4
#define PGPLOT_CT_AIPS 5

/* the different Stokes parameters we can make */
#define STOKES_I -1
#define STOKES_Q -2
#define STOKES_U -3
#define STOKES_V -4

/* some PGPLOT constants that I should really put into a global header
   file for future use */
#define PGSCI_BACKGROUND 0
#define PGSCI_INVERSE 1
#define PGSCI_RED 2
#define PGSCI_GREEN 3
#define PGSCI_BLUE 4
#define PGSCI_CYAN 5
#define PGSCI_MAGENTA 6
#define PGSCI_YELLOW 7
#define PGSCI_ORANGE 8
#define PGSCI_GREENYELLOW 9
#define PGSCI_GREENCYAN 10
#define PGSCI_BLUECYAN 11
#define PGSCI_BLUEMAGENTA 12
#define PGSCI_REDMAGENTA 13
#define PGSCI_DARKGREY 14
#define PGSCI_LIGHTGREY 15

#define PGSLS_FULL 1
#define PGSLS_DASHED 2
#define PGSLS_DOTDASH 3
#define PGSLS_DOTTED 4
#define PGSLS_DASHDOTS 5

#define PGLEN_NORMALIZED_DEVICE 0
#define PGLEN_INCHES 1
#define PGLEN_MM 2
#define PGLEN_ABSOLUTE_DEVICE 3
#define PGLEN_WORLD 4
#define PGLEN_VIEWPORT_FRACTION 5

#define PGQINF_ITEM_VERSION "VERSION"
#define PGQINF_ITEM_STATE "STATE"
#define PGQINF_ITEM_USER "USER"
#define PGQINF_ITEM_DATETIME "NOW"
#define PGQINF_ITEM_DEVICE "DEVICE"
#define PGQINF_ITEM_TYPE "TYPE"
#define PGQINF_ITEM_DEVICETYPE "DEV/TYPE"
#define PGQINF_ITEM_ISHARDCOPY "HARDCOPY"
#define PGQINF_ITEM_ISTERMINAL "TERMINAL"
#define PGQINF_ITEM_HASCURSOR "CURSOR"
#define PGQINF_ITEM_HASSCROLL "SCROLL"
#define PGQINF_VALUE_UNSUCCESFUL "?"
#define PGQINF_VALUE_YES "YES"
#define PGQINF_VALUE_NO "NO"
#define PGQINF_VALUE_OPEN "OPEN"
#define PGQINF_VALUE_CLOSED "CLOSED"

/* how should the peak summary plot be divided up?
   the following settings are to make a 5x6 array of plots */
#define PANEL_MIN_ANT 1
#define PANEL_MAX_ANT 6
#define PANEL_X_COLS 5
#define PANEL_Y_ROWS 6
#define PANEL_MIN_IF 1
#define PANEL_MAX_IF 2

#define PGPLOT_MEDIUM_SCREEN -1
#define PGPLOT_MEDIUM_PAPER -2

/* make a safe malloc, realloc, calloc & free */
#define MALLOC(a,b) if ((a=malloc((b)))==NULL) { fprintf(stderr,"malloc failed at %s:%d!\n",__FILE__,__LINE__); exit(-1); }
#define REALLOC(a,b) if ((a=realloc(a,(b)))==NULL) { fprintf(stderr,"realloc failed at %s:%d!\n",__FILE__,__LINE__); exit(-1); }
#define CALLOC(a,b,c) if ((a=calloc(b,c))==NULL) { fprintf(stderr,"calloc failed at %s:%d!\n",__FILE__,__LINE__); exit(-1); }
#define FREE(a) { free(a); a=NULL; }

/* matching criteria are set here and checked programmatically 
   making it easier to tweak what we're looking for */
/* the number of instruments that can be considered */
#define CONSIDER_N_INSTRUMENTS 1
/* the list of instruments that can be considered */
char consider_instrument[CONSIDER_N_INSTRUMENTS][BUFSIZE]={"ATCA"};

/* the minimum number of IFs required */
#define MIN_N_IFS 2
/* uncomment this next line if this is a minimum number instead of an exact
   number */
/* #define MORE_IFS */

/* uncomment this next line if you don't care about simultaneity */
/* #define NOT_SIMULTANEOUS */

/* the maximum number of frequency channels */
#define MAX_N_CHANNELS 32
/* uncomment this next line if this is a maximum number instead of an exact
   number */
/* #define LESS_CHANNELS */

/* the maximum frequency to consider, in Hz */
#define MAX_IF_FREQUENCY 1600E6

/* the minimum bandwidth, in Hz */
#define MIN_BANDWIDTH 64E6

/* don't count visibilities within this number of from the scan boundaries */
#define SCAN_BOUNDARY_LIMIT 5
/* uncomment this line if you want to keep data at the end of scans */
#define KEEP_ENDOFSCAN

typedef struct xmlfile_info_ {
  char program_date[BUFSIZE];
  char program_time[BUFSIZE];
  char timezone[BUFSIZE];
  int nearby_spectra;
  float minimum_significance;
  float baseline_match_fraction;
} xmlfile_info;

typedef struct peak_baseline_ {
  int if_number;
  int baseline_ant1;
  int baseline_ant2;
  float baseline_length;
  int lag_channel;
  float baseline_lag;
  float baseline_lag_err;
  float amplitude;
  float spectrum_average_amplitude;
  float spectrum_rms_noise;
  float significance;
} peak_baseline;

typedef struct transient_peak_ {
  int peak_number;
  char source_name[BUFSIZE];
  char string_observation_time[BUFSIZE];
  time_t observation_time;
  int n_lags;
  peak_baseline **lag_peaks;
  int available_baselines;
  int peak_baselines;
  float peak_baselines_fraction;
  int transient_candidate;
} transient_peak;

typedef struct archive_file_header_ {
  int header_number;
  char obs_date[BUFSIZE];
  int cycle_time;
  int number_ifs;
  int number_sources;
  int *if_number;
  float *if_central_frequency;
  float *if_bandwidth;
  int *if_channels;
  int *source_number;
  char **source_name;
  int usable;
  float usable_time;
} archive_file_header;

typedef struct archive_file_ {
  char path[BUFSIZE];
  char filename[BUFSIZE];
  int usable;
  float usable_time;
  int n_peaks;
  transient_peak **transient_peaks;
  int n_plots;
  char **plot_peak_number;
  char **plot_name;
  int n_headers;
  archive_file_header **headers;
  char benchmark_file_read_time[BUFSIZE];
  char benchmark_fourier_transform_time[BUFSIZE];
  char benchmark_polarisation_combination_time[BUFSIZE];
  char benchmark_peak_search_time[BUFSIZE];
  char benchmark_peak_verification_time[BUFSIZE];
  struct archive_file_ *next;
} archive_file;

typedef struct spectral_info_ {
  float average_amplitude;
  float rms_noise;
  int npeaks_found;
  int *peak_channels;
  int *peak_verified;
  int *peak_number;
  float estimated_average_amplitude;
  float estimated_average_amplitude_error;
  float estimated_rms_noise;
  double fit_c0;
  double fit_c1;
} spectral_info;

/* the structure to hold all the data in */
typedef struct correlation_data_ {
  float *vis;
  float *wgt;
  int baseline;
  int islongest_baseline;
  float ut;
  time_t observation_time;
  float u;
  float v;
  float w;
  int flag;
  int bin;
  int if_no;
  int sourceno;
  char sourcename[BUFSIZE];
  char **stokes_name;
  int if_nfreq;
  int if_nstok;
  float if_bandwidth;
  int scannum;
  float *pols_combined;
  float **individual_pols;
  char combined_pol_name[BUFSIZE];
  spectral_info spectrum_info;
  struct correlation_data_ *next;
  struct correlation_data_ *prev;
  struct correlation_data_ *separate_next;
  struct correlation_data_ *samebaseline_next;
  struct correlation_data_ *samebaseline_prev;
  struct correlation_data_ *raw_data;
} correlation_data;

void strip_trailing_spaces(char *string);
archive_file *file_info(char *filename,archive_file **file_list);
void write_xml_file(char *report_file,archive_file **file_list);
void build_file_database(char *report_file,archive_file **file_list);
void free_archive_file (archive_file **unwanted_file);
int peak_summary_plot(correlation_data **all_correlations,int interactive,char *pgplot_device);
int want_this_file(int *longest_baseline);
int clear_correlations(correlation_data **all_correlations);
int remove_correlation(correlation_data **correlation_to_remove);
int transient_peak_menu(archive_file **current_file,int current_peak,int interactive);
int read_data(correlation_data **all_correlations,int longest_baseline,int *keep_reading);
int menu_question(char *question,int n_options,int default_option,...);
int change_to_lagspace(correlation_data **all_correlations,
		       correlation_data **all_transforms);
void copy_correlation_info(correlation_data **original,
			   correlation_data **destination);
int combine_polarisations(correlation_data **all_correlations,int stokes_to_make);
void which_panel(int baseline,int if_no,int *xpanel,int *ypanel);
int add_correlation(correlation_data **correlation_list,float **vis,float **wgt,int baseline,
		    float ut,float u,float v,float w,int flag,int bin,int if_no,int sourceno,
		    int if_nfreq,int if_nstok,int longest_baseline,int scannum);
float standard_deviation(float **values,int nelements,float *average);
float meanvalue(float **values,int nelements);
void baseline_to_antennae(int baseline_number,int *antenna1,int *antenna2);
void strip_path(char *full_path,char *filename);
xmlChar *ConvertInput(const char *in, const char *encoding);

char currently_loaded_file[BUFSIZE];
correlation_data *current_file_correlations;
xmlfile_info info_xmlfile;

/* now the argp command-line argument parsing stuff */
static struct argp_option options[] = {
  {"xml-file",'f',"XML_FILE",0,"the name of an XML report file from transient_search_data"},
  {"make-file-list",'m',0,0,"make a file summarising all the files in the xml-file"},
  {"search-for-file",'s',"SEARCH_FILE",0,"search for a file and summarise its properties"},
  {"inspect-transients",'i',0,0,"interactive inspection mode of transients"},
  {"output-filename",'o',"OUTPUT_FILENAME",0,"the name of the file to output"},
  {"useful-files-only",'u',0,0,"include only files with useful data in the file list output"},
  {"make-plots",'p',0,0,"make plots of each potential transient"},
  {"individual-limits",'l',0,0,"use the best limits for each individual plot panel"},
  { 0 }
};

struct arguments {
  int n_xml_files;
  char **xml_files;
  int inspect_transients;
  int make_file_list;
  int useful_files_only;
  int make_plots;
  int search_for_file;
  int n_files_to_find;
  char **files_to_find;
  char output_filename[BUFSIZE];
  int individual_limits;
};

struct arguments arguments;

static error_t parse_opt(int key,char *arg,struct argp_state *state){
  struct arguments *arguments = state->input;

  switch(key){
  case 'f': /* an XML report file */
    arguments->n_xml_files++;
    REALLOC(arguments->xml_files,arguments->n_xml_files*sizeof(char*));
    MALLOC(arguments->xml_files[arguments->n_xml_files-1],BUFSIZE*sizeof(char));
    strncpy(arguments->xml_files[arguments->n_xml_files-1],arg,BUFSIZE);
    break;
  case 'm': /* make a file listing the files in the XML report files */
    arguments->make_file_list=YES;
    break;
  case 's': /* search for some file */
    arguments->search_for_file=YES;
    arguments->n_files_to_find++;
    REALLOC(arguments->files_to_find,arguments->n_files_to_find*sizeof(char*));
    MALLOC(arguments->files_to_find[arguments->n_files_to_find-1],
	   BUFSIZE*sizeof(char));
    strncpy(arguments->files_to_find[arguments->n_files_to_find-1],arg,BUFSIZE);
    break;
  case 'i': /* inspect the detected transient candidates */
    arguments->inspect_transients=YES;
    break;
  case 'o': /* the desired name of the output file */
    strncpy(arguments->output_filename,arg,BUFSIZE);
    break;
  case 'u': /* make a list of useful files only */
    arguments->useful_files_only=YES;
    break;
  case 'p': /* make the transient peak plots */
    arguments->make_plots=YES;
    break;
  case 'l': /* use individual plot limits */
    arguments->individual_limits=YES;
    break;
  }

  return(0);
}

static struct argp argp = { options, parse_opt, 0, 0 };

void strip_trailing_spaces(char *string){
  /* this routine takes a string that may or may not have trailing spaces before the null 
     terminator, and removes those trailing spaces if they exist */
  int i;
  
  for (i=strlen(string)-1;i>=0;i--){
    if (string[i]==' '){
      string[i]='\0'; /* replace the space with the null terminator */
    } else {
      /* not a space so we stop here */
      break;
    }
  }
}

archive_file *file_info(char *filename,archive_file **file_list){
  archive_file *cycle_file=NULL;

  for (cycle_file=*file_list;cycle_file!=NULL;
       cycle_file=cycle_file->next){
    if (strstr(filename,cycle_file->filename)!=NULL)
      return(cycle_file);
  }
  return(NULL);
}

#define XMLREPORT(routine,...) rc=routine(__VA_ARGS__);			\
  if (rc<0) { printf("XML writing error at %s:%d\n",__FILE__,__LINE__);	}
#define XMLCONVERT(string,xmlstring,...) if (xmlstring!=NULL){		\
    xmlFree(xmlstring); } snprintf(string,BUFSIZE,__VA_ARGS__);		\
  xmlstring=ConvertInput(string,MY_ENCODING);

void write_xml_file(char *report_file,archive_file **file_list){
  xmlTextWriterPtr xml_report_writer;
  int rc,i,j;
  char tmp[BUFSIZE];
  xmlChar *xmltmp=NULL;
  archive_file *cycle_file=NULL;

  xml_report_writer=xmlNewTextWriterFilename(report_file,0);
  if (xml_report_writer==NULL){
    printf("unable to open %s for writing!\n",report_file);
    return;
  }

  /* start the XML report with the XML default for the version, encoding and standalone
     declaration */
  XMLREPORT(xmlTextWriterStartDocument,xml_report_writer,NULL,MY_ENCODING,NULL);
  /* start the root element, call it ATCAARCHIVE */
  XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "ATCAARCHIVE");
  /* start an element called DATETIME as a child of ATCAARCHIVE */
  XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "DATETIME");
  /* now give it the attribute date and value */
  XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,BAD_CAST "date",
	    BAD_CAST info_xmlfile.program_date);
  XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,BAD_CAST "time",
	    BAD_CAST info_xmlfile.program_time);
  XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,BAD_CAST "timezone",
	    BAD_CAST info_xmlfile.timezone);
  /* comments must come after attributes */
  XMLCONVERT(tmp,xmltmp,"the date and time at program start");
  XMLREPORT(xmlTextWriterWriteComment,xml_report_writer,xmltmp);
  /* close the element DATETIME */
  XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
  /* now write out the parameters we're using */
  XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "PARAMETERS");
  snprintf(tmp,BUFSIZE,"%d",info_xmlfile.nearby_spectra);
  XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	    BAD_CAST "nearby_spectra",BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%.1f",info_xmlfile.minimum_significance);
  XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	    BAD_CAST "minimum_significance", BAD_CAST tmp);
  snprintf(tmp,BUFSIZE,"%.3f",info_xmlfile.baseline_match_fraction);
  XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	    BAD_CAST "baseline_match_fraction",BAD_CAST tmp);
  XMLCONVERT(tmp,xmltmp,"nearby_spectra is the number of nearby spectra use to estimate the average power and RMS noise of a target spectrum");
  XMLREPORT(xmlTextWriterWriteComment,xml_report_writer,xmltmp);
  XMLCONVERT(tmp,xmltmp,"minimum_significance is the minimum number of RMS noise levels above the average a peak needs to be to be considered");
  XMLREPORT(xmlTextWriterWriteComment,xml_report_writer,xmltmp);
  XMLCONVERT(tmp,xmltmp,"baseline_match_fraction is the fraction of baselines that must have a corresponding peak to be considered");
  XMLREPORT(xmlTextWriterWriteComment,xml_report_writer,xmltmp);
  XMLREPORT(xmlTextWriterEndElement,xml_report_writer);

  /* now loop over the files */
  for (cycle_file=*file_list;cycle_file!=NULL;cycle_file=cycle_file->next){
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "FILE");
    snprintf(tmp,BUFSIZE,"%s/%s",cycle_file->path,cycle_file->filename);
    XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	      BAD_CAST "filename",BAD_CAST tmp);
    /* now loop over the headers in the file */
    for (i=0;i<cycle_file->n_headers;i++){
      XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "HEADER");
      snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->header_number);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "number",BAD_CAST tmp);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "date",BAD_CAST cycle_file->headers[i]->obs_date);
      snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->cycle_time);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "cycle_time",BAD_CAST tmp);
      snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->number_ifs);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "number_ifs",BAD_CAST tmp);
      snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->number_sources);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "number_sources",BAD_CAST tmp);
      /* loop over the IFs */
      for (j=0;j<cycle_file->headers[i]->number_ifs;j++){
	if (j>0){
	  /* have to close the last IF element first */
	  XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
	}
	XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "IF");
	snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->if_number[j]);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "number",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.1f",cycle_file->headers[i]->if_central_frequency[j]);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "central_frequency",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.1f",cycle_file->headers[i]->if_bandwidth[j]);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "bandwidth",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->if_channels[j]);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "channels",BAD_CAST tmp);
      }
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
      /* loop over the sources */
      for (j=0;j<cycle_file->headers[i]->number_sources;j++){
	if (j>0){
	  /* have to close the last IF element first */
	  XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
	}
	XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "SOURCE");
	snprintf(tmp,BUFSIZE,"%d",cycle_file->headers[i]->source_number[j]);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "number",BAD_CAST tmp);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "name",BAD_CAST cycle_file->headers[i]->source_name[j]);
      }
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
      /* is this data usable */
      if (cycle_file->headers[i]->usable==YES){
	XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "SUMMARY");
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "usable",BAD_CAST "yes");
	snprintf(tmp,BUFSIZE,"%.1f",cycle_file->headers[i]->usable_time);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "usable_time",BAD_CAST tmp);
	XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
      } else if (cycle_file->headers[i]->usable==NO){
	XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "SUMMARY");
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "usable",BAD_CAST "no");
	XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
      }
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    }
    /* now write out the file reading time benchmark */
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	      BAD_CAST "file_read_time",
	      BAD_CAST cycle_file->benchmark_file_read_time);
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    /* now a file usability summary */
    if (cycle_file->usable==YES){
      XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "SUMMARY");
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "usable",BAD_CAST "yes");
      snprintf(tmp,BUFSIZE,"%.1f",cycle_file->usable_time);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "total_usable_time",BAD_CAST tmp);
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    } else if (cycle_file->usable==NO){
      XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "SUMMARY");
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "usable",BAD_CAST "no");
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    }
    /* now the other benchmarks */
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	      BAD_CAST "fourier_transform_time",
	      BAD_CAST cycle_file->benchmark_fourier_transform_time);
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	      BAD_CAST "polarisation_combination_time",
	      BAD_CAST cycle_file->benchmark_polarisation_combination_time);
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	      BAD_CAST "peak_search_time",
	      BAD_CAST cycle_file->benchmark_peak_search_time);
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    /* now go through all the peaks that were checked */
    for (i=0;i<cycle_file->n_peaks;i++){
      XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "PEAK");
      snprintf(tmp,BUFSIZE,"%04d",cycle_file->transient_peaks[i]->peak_number);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "number",BAD_CAST tmp);
      /* now go through all the lags attributed to this peak */
      for (j=0;j<cycle_file->transient_peaks[i]->n_lags;j++){
	XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "FOUND");
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "source_name",
		  BAD_CAST cycle_file->transient_peaks[i]->source_name);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "observation_time",
		  BAD_CAST cycle_file->transient_peaks[i]->string_observation_time);
	snprintf(tmp,BUFSIZE,"%d",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->if_number);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "if_number",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%d",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_ant1);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "baseline_ant1",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%d",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_ant2);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "baseline_ant2",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_length);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "baseline_length",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%d",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->lag_channel);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "lag_channel",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_lag);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "baseline_lag",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_lag_err);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "baseline_lag_err",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->amplitude);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "amplitude",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->spectrum_average_amplitude);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "spectrum_average_amplitude",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->spectrum_rms_noise);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "spectrum_rms_noise",BAD_CAST tmp);
	snprintf(tmp,BUFSIZE,"%.2f",
		 cycle_file->transient_peaks[i]->lag_peaks[j]->significance);
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "significance",BAD_CAST tmp);
	XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
      }
      /* now summarise the peaks/baseline situation */
      XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "SUMMARY");
      snprintf(tmp,BUFSIZE,"%d",
	       cycle_file->transient_peaks[i]->available_baselines);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "available_baselines",BAD_CAST tmp);
      snprintf(tmp,BUFSIZE,"%d",
	       cycle_file->transient_peaks[i]->peak_baselines);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "peak_baselines",BAD_CAST tmp);
      snprintf(tmp,BUFSIZE,"%.4f",
	       cycle_file->transient_peaks[i]->peak_baselines_fraction);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "peak_baselines_fraction",BAD_CAST tmp);
      if (cycle_file->transient_peaks[i]->transient_candidate==YES){
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "transient_candidate",BAD_CAST "yes");
      } else if (cycle_file->transient_peaks[i]->transient_candidate==NO){
	XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		  BAD_CAST "transient_candidate",BAD_CAST "no");
      }
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    }
    /* another benchmark */
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "BENCHMARK");
    XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
	      BAD_CAST "peak_verification_time",
	      BAD_CAST cycle_file->benchmark_peak_verification_time);
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    /* now the plot info */
    XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "PLOT");
    for (i=0;i<cycle_file->n_plots;i++){
      XMLREPORT(xmlTextWriterStartElement,xml_report_writer,BAD_CAST "PEAK");
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "number",
		BAD_CAST cycle_file->plot_peak_number[i]);
      XMLREPORT(xmlTextWriterWriteAttribute,xml_report_writer,
		BAD_CAST "filename",
		BAD_CAST cycle_file->plot_peak_number[i]);
      XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    }
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
    XMLREPORT(xmlTextWriterEndElement,xml_report_writer);
  }
  XMLREPORT(xmlTextWriterEndDocument,xml_report_writer);
  xmlFreeTextWriter(xml_report_writer);
}

void build_file_database(char *report_file,archive_file **file_list){
  xmlTextReaderPtr reader;
  xmlChar *name=NULL,*attrname=NULL,*attrvalue=NULL;
  int ret,nodetype,slash_offset,depth,i,current_if_number,current_source_number;
  int in_plot=NO;
  archive_file *new_file=NULL;
  transient_peak *new_peak=NULL;
  struct tm observation_time;
  peak_baseline *new_baseline=NULL;
  archive_file_header *new_header=NULL;

  reader=xmlReaderForFile(report_file,"ISO-8859-1",
			  XML_PARSE_RECOVER|XML_PARSE_NOENT|XML_PARSE_XINCLUDE);
  if (reader==NULL){
    printf(" %s: cannot open file\n",report_file);
    return;
  }

  while((ret=xmlTextReaderRead(reader))==1){
    name=xmlTextReaderName(reader);
    nodetype=xmlTextReaderNodeType(reader);
    depth=xmlTextReaderDepth(reader);
    if (strcmp((char*)name,"FILE")==0){
      if (nodetype==XML_READER_TYPE_ELEMENT){
	/* the start of a new file */
	/* make a new file structure */
	new_file=malloc(sizeof(archive_file));
	new_file->transient_peaks=NULL;
	new_file->next=NULL;
	new_file->n_peaks=0;
	new_file->n_plots=0;
	new_file->plot_peak_number=NULL;
	new_file->plot_name=NULL;
	new_file->n_headers=0;
	new_file->headers=NULL;
	/* get the name of the file from the attributes */
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  if (strcmp((char*)attrname,"filename")==0){
	    attrvalue=xmlTextReaderValue(reader);
	    if (attrvalue!=NULL){
	      for (slash_offset=strlen((char*)attrvalue);
		   slash_offset>0;slash_offset--){
		if (attrvalue[slash_offset]=='/'){
		  break;
		}
	      }
	      if (new_file!=NULL){
		strncpy(new_file->path,(char*)attrvalue,slash_offset);
		if (slash_offset!=0)
		  strcpy(new_file->filename,(char*)attrvalue+slash_offset+1);
		else
		  strcpy(new_file->filename,(char*)attrvalue);
	      }
	      xmlFree(attrvalue);
	    }
	  }
	  xmlFree(attrname);
	}
      } else if (nodetype==XML_READER_TYPE_END_ELEMENT){
	/* the end of a file */
	/* add the file to the head of the list */
	if (new_file!=NULL){
	  new_file->next=*file_list;
	  *file_list=new_file;
	}
	new_file=NULL;
      }
    } else if (strcmp((char*)name,"SUMMARY")==0){
      if (depth==2){
	/* a summary of if there is any usable data in the file */
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  if (strcmp((char*)attrname,"usable")==0){
	    attrvalue=xmlTextReaderValue(reader);
	    if (attrvalue!=NULL){
	      if (strcmp((char*)attrvalue,"yes")==0){
		if (new_file!=NULL)
		  new_file->usable=YES;
	      } else if (strcmp((char*)attrvalue,"no")==0){
		if (new_file!=NULL)
		  new_file->usable=NO;
	      }
	      xmlFree(attrvalue);
	    }
	  } else if (strcmp((char*)attrname,"total_usable_time")==0){
	    attrvalue=xmlTextReaderValue(reader);
	    new_file->usable_time=atof((char*)attrvalue);
	    xmlFree(attrvalue);
	  }
	  xmlFree(attrname);
	}
      } else if (depth==3){
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  if (new_peak!=NULL){
	    attrname=xmlTextReaderName(reader);
	    attrvalue=xmlTextReaderValue(reader);
	    if (strcmp((char*)attrname,"available_baselines")==0){
	      new_peak->available_baselines=atoi((char*)attrvalue);
	    } else if (strcmp((char*)attrname,"peak_baselines")==0){
	      new_peak->peak_baselines=atoi((char*)attrvalue);
	    } else if (strcmp((char*)attrname,"peak_baselines_fraction")==0){
	      new_peak->peak_baselines_fraction=atof((char*)attrvalue);
	    } else if (strcmp((char*)attrname,"transient_candidate")==0){
	      if (strcmp((char*)attrvalue,"yes")==0){
		new_peak->transient_candidate=YES;
	      } else if (strcmp((char*)attrvalue,"no")==0){
		new_peak->transient_candidate=NO;
	      }
	    }
	    xmlFree(attrname);
	    xmlFree(attrvalue);
	  } else if (new_header!=NULL){
	    /* we find out whether the header describes usable data */
	    attrname=xmlTextReaderName(reader);
	    attrvalue=xmlTextReaderValue(reader);
	    if (strcmp((char*)attrname,"usable")==0){
	      if (strcmp((char*)attrvalue,"yes")==0)
		new_header->usable=YES;
	      else if (strcmp((char*)attrvalue,"no")==0)
		new_header->usable=NO;
	    } else if (strcmp((char*)attrname,"usable_time")==0){
	      new_header->usable_time=atof((char*)attrvalue);
	    }
	    xmlFree(attrname);
	    xmlFree(attrvalue);
	  }
	}
      }
    } else if (strcmp((char*)name,"PEAK")==0){
      if (in_plot==NO){
	if (nodetype==XML_READER_TYPE_ELEMENT){
	  /* a peak number */
	  new_peak=malloc(sizeof(transient_peak));
	  new_peak->n_lags=0;
	  new_peak->lag_peaks=NULL;
	  while(xmlTextReaderMoveToNextAttribute(reader)==1){
	    attrname=xmlTextReaderName(reader);
	    attrvalue=xmlTextReaderValue(reader);
	    if (attrvalue==NULL){
	      /* oops this is a bad thing! */
	      xmlFree(attrname);
	      free(new_peak);
	      new_peak=NULL;
	      break;
	    }
	    if (strcmp((char*)attrname,"number")==0){
	      new_peak->peak_number=atoi((char*)attrvalue);
	    }
	    xmlFree(attrname);
	    xmlFree(attrvalue);
	  }
	} else if (nodetype==XML_READER_TYPE_END_ELEMENT){
	  /* attach the peak to the file structure */
	  if (new_peak!=NULL){
	    new_file->n_peaks++;
	    new_file->transient_peaks=realloc(new_file->transient_peaks,
					      new_file->n_peaks*sizeof(transient_peak));
	    new_file->transient_peaks[new_file->n_peaks-1]=new_peak;
	    new_peak=NULL;
	  }
	}
      } else if (in_plot==YES){
	/* this is telling us about the file the peak number is 
	   plotted in */
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (strcmp((char*)attrname,"number")==0){
	    new_file->n_plots++;
	    REALLOC(new_file->plot_peak_number,
		    new_file->n_plots*sizeof(char*));
	    REALLOC(new_file->plot_name,
		    new_file->n_plots*sizeof(char*));
	    MALLOC(new_file->plot_peak_number[new_file->n_plots-1],
		   BUFSIZE*sizeof(char));
	    MALLOC(new_file->plot_name[new_file->n_plots-1],
		   BUFSIZE*sizeof(char));
	    strcpy(new_file->plot_peak_number[new_file->n_plots-1],
		   (char*)attrvalue);
	  } else if (strcmp((char*)attrname,"filename")==0){
	    strcpy(new_file->plot_name[new_file->n_plots-1],(char*)attrvalue);
	  }
	  xmlFree(attrname);
	  xmlFree(attrvalue);
	}
      }
    } else if (strcmp((char*)name,"FOUND")==0){
      /* we've got a peak */
      if (new_peak!=NULL){
	new_baseline=malloc(sizeof(peak_baseline));
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (attrvalue==NULL){
	    /* oops this is a bad thing! */
	    xmlFree(attrname);
	    free(new_peak);
	    new_peak=NULL;
	    break;
	  }
	  if (strcmp((char*)attrname,"source_name")==0){
	    strcpy(new_peak->source_name,(char*)attrvalue);
	  } else if (strcmp((char*)attrname,"observation_time")==0){
	    strcpy(new_peak->string_observation_time,(char*)attrvalue);
	    sscanf((char*)attrvalue,"%2d:%2d:%2d %2d/%2d/%4d",
		   &(observation_time.tm_hour),&(observation_time.tm_min),
		   &(observation_time.tm_sec),&(observation_time.tm_mday),
		   &(observation_time.tm_mon),&(observation_time.tm_year));
	    observation_time.tm_mon--;
	    observation_time.tm_year-=1900;
	    new_peak->observation_time=timegm(&observation_time);
	  } else if (strcmp((char*)attrname,"if_number")==0){
	    new_baseline->if_number=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_ant1")==0){
	    new_baseline->baseline_ant1=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_ant2")==0){
	    new_baseline->baseline_ant2=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_length")==0){
	    new_baseline->baseline_length=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"lag_channel")==0){
	    new_baseline->lag_channel=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_lag")==0){
	    new_baseline->baseline_lag=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_lag_err")==0){
	    new_baseline->baseline_lag_err=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"amplitude")==0){
	    new_baseline->amplitude=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"spectrum_average_amplitude")==0){
	    new_baseline->spectrum_average_amplitude=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"spectrum_rms_noise")==0){
	    new_baseline->spectrum_rms_noise=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"significance")==0){
	    new_baseline->significance=atof((char*)attrvalue);
	    new_peak->n_lags++;
	    new_peak->lag_peaks=realloc(new_peak->lag_peaks,
					new_peak->n_lags*sizeof(peak_baseline));
	    new_peak->lag_peaks[new_peak->n_lags-1]=new_baseline;
	    new_baseline=NULL;
	  }
	  xmlFree(attrname);
	  xmlFree(attrvalue);
	}
      }
    } else if (strcmp((char*)name,"DATETIME")==0){
      if (nodetype==XML_READER_TYPE_ELEMENT){
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (strcmp((char*)attrname,"date")==0){
	    strcpy(info_xmlfile.program_date,(char*)attrvalue);
	  } else if (strcmp((char*)attrname,"time")==0){
	    strcpy(info_xmlfile.program_time,(char*)attrvalue);
	  } else if (strcmp((char*)attrname,"timezone")==0){
	    strcpy(info_xmlfile.timezone,(char*)attrvalue);
	  }
	  xmlFree(attrname);
	  xmlFree(attrvalue);
	}
      }
    } else if (strcmp((char*)name,"PARAMETERS")==0){
      if (nodetype==XML_READER_TYPE_ELEMENT){
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (strcmp((char*)attrname,"nearby_spectra")==0){
	    info_xmlfile.nearby_spectra=atoi((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"minimum_significance")==0){
	    info_xmlfile.minimum_significance=atof((char*)attrvalue);
	  } else if (strcmp((char*)attrname,"baseline_match_fraction")==0){
	    info_xmlfile.baseline_match_fraction=atof((char*)attrvalue);
	  }
	  xmlFree(attrname);
	  xmlFree(attrvalue);
	}
      }
    } else if (strcmp((char*)name,"HEADER")==0){
      if (nodetype==XML_READER_TYPE_ELEMENT){
	/* allocate a new header */
	MALLOC(new_header,sizeof(archive_file_header));
	while(xmlTextReaderMoveToNextAttribute(reader)==1){
	  attrname=xmlTextReaderName(reader);
	  attrvalue=xmlTextReaderValue(reader);
	  if (strcmp((char*)attrname,"number")==0)
	    new_header->header_number=atoi((char*)attrvalue);
	  else if (strcmp((char*)attrname,"date")==0)
	    strcpy(new_header->obs_date,(char*)attrvalue);
	  else if (strcmp((char*)attrname,"cycle_time")==0)
	    new_header->cycle_time=atoi((char*)attrvalue);
	  else if (strcmp((char*)attrname,"number_ifs")==0){
	    new_header->number_ifs=atoi((char*)attrvalue);
	    MALLOC(new_header->if_number,new_header->number_ifs*sizeof(int));
	    MALLOC(new_header->if_central_frequency,
		   new_header->number_ifs*sizeof(float));
	    MALLOC(new_header->if_bandwidth,
		   new_header->number_ifs*sizeof(float));
	    MALLOC(new_header->if_channels,
		   new_header->number_ifs*sizeof(int));
	  } else if (strcmp((char*)attrname,"number_sources")==0){
	    new_header->number_sources=atoi((char*)attrvalue);
	    MALLOC(new_header->source_number,
		   new_header->number_sources*sizeof(int));
	    MALLOC(new_header->source_name,
		   new_header->number_sources*sizeof(char*));
	    for (i=0;i<new_header->number_sources;i++)
	      MALLOC(new_header->source_name[i],BUFSIZE*sizeof(char));
	  }
	  xmlFree(attrname);
	  xmlFree(attrvalue);
	}
	current_if_number=0;
	current_source_number=0;
      } else if (nodetype==XML_READER_TYPE_END_ELEMENT){
	/* we've finished reading about the header, so we add this
	   info to the file structure */
	new_file->n_headers++;
	REALLOC(new_file->headers,
		new_file->n_headers*sizeof(archive_file_header*));
	new_file->headers[new_file->n_headers-1]=new_header;
	new_header=NULL;
      }
    } else if (strcmp((char*)name,"IF")==0){
      /* we're in a header element, so we continue filling in the 
         new header structure */
      while(xmlTextReaderMoveToNextAttribute(reader)==1){
	attrname=xmlTextReaderName(reader);
	attrvalue=xmlTextReaderValue(reader);
	if (strcmp((char*)attrname,"number")==0)
	  new_header->if_number[current_if_number]=atoi((char*)attrvalue);
	else if (strcmp((char*)attrname,"central_frequency")==0)
	  new_header->if_central_frequency[current_if_number]=
	    atof((char*)attrvalue);
	else if (strcmp((char*)attrname,"bandwidth")==0)
	  new_header->if_bandwidth[current_if_number]=
	    atof((char*)attrvalue);
	else if (strcmp((char*)attrname,"channels")==0)
	  new_header->if_channels[current_if_number]=
	    atoi((char*)attrvalue);
	xmlFree(attrname);
	xmlFree(attrvalue);
      }
      current_if_number++;
    } else if (strcmp((char*)name,"SOURCE")==0){
      /* we're in a header element, so we continue filling in the 
         new header structure */
      while(xmlTextReaderMoveToNextAttribute(reader)==1){
	attrname=xmlTextReaderName(reader);
	attrvalue=xmlTextReaderValue(reader);
	if (strcmp((char*)attrname,"number")==0)
	  new_header->source_number[current_source_number]=
	    atoi((char*)attrvalue);
	else if (strcmp((char*)attrname,"name")==0)
	  strcpy(new_header->source_name[current_source_number],
		 (char*)attrvalue);
	xmlFree(attrname);
	xmlFree(attrvalue);
      }
      current_source_number++;
    } else if (strcmp((char*)name,"BENCHMARK")==0){
      /* read in how fast the code did all the stuff */
      while(xmlTextReaderMoveToNextAttribute(reader)==1){
	attrname=xmlTextReaderName(reader);
	attrvalue=xmlTextReaderValue(reader);
	if (strcmp((char*)attrname,"file_read_time")==0)
	  strcpy(new_file->benchmark_file_read_time,
		 (char*)attrvalue);
	else if (strcmp((char*)attrname,"fourier_transform_time")==0)
	  strcpy(new_file->benchmark_fourier_transform_time,
		 (char*)attrvalue);
	else if (strcmp((char*)attrname,"polarisation_combination_time")==0)
	  strcpy(new_file->benchmark_polarisation_combination_time,
		 (char*)attrvalue);
	else if (strcmp((char*)attrname,"peak_search_time")==0)
	  strcpy(new_file->benchmark_peak_search_time,
		 (char*)attrvalue);
	else if (strcmp((char*)attrname,"peak_verification_time")==0)
	  strcpy(new_file->benchmark_peak_verification_time,
		 (char*)attrvalue);
	xmlFree(attrname);
	xmlFree(attrvalue);
      }
    } else if (strcmp((char*)name,"PLOT")==0){
      /* details about the plot the program made */
      if (nodetype==XML_READER_TYPE_ELEMENT){
	in_plot=YES;
      } else if (nodetype==XML_READER_TYPE_END_ELEMENT){
	in_plot=NO;
      }
    }
  }
  xmlFreeTextReader(reader);
}

void free_archive_file (archive_file **unwanted_file) {
  int i,j;
  for (i=0;i<(*unwanted_file)->n_peaks;i++){
    for (j=0;j<(*unwanted_file)->transient_peaks[i]->n_lags;j++){
      free((*unwanted_file)->transient_peaks[i]->lag_peaks[j]);
    }
    free((*unwanted_file)->transient_peaks[i]->lag_peaks);
  }
  free((*unwanted_file)->transient_peaks);
  free(*unwanted_file);
}

int main(int argc,char *argv[]){
  char filename[BUFSIZE],timestring[BUFSIZE],autosave_filename[BUFSIZE];
  int transient_count,i,j,nf=0,option_made,xml_changed=NO;
  archive_file *file_list=NULL,*cycle_file=NULL,*free_file=NULL;
  FILE *list_of_files=NULL;

  /* initialise the arguments */
  arguments.n_xml_files=0;
  arguments.xml_files=NULL;
  arguments.inspect_transients=NO;
  arguments.make_file_list=NO;
  arguments.make_plots=NO;
  arguments.useful_files_only=NO;
  arguments.search_for_file=NO;
  arguments.n_files_to_find=0;
  arguments.files_to_find=NULL;
  arguments.output_filename[0]='\0';
  arguments.individual_limits=NO;

  /* parse the command line arguments */
  argp_parse(&argp,argc,argv,0,0,&arguments);

  if ((arguments.inspect_transients==NO)&&
      (arguments.make_file_list==NO)&&
      (arguments.search_for_file==NO)&&
      (arguments.make_plots==NO)){
    /* no mode has been selected, so we default to inspecting
       the detected transients */
    arguments.inspect_transients=YES;
  }

  current_file_correlations=NULL; /* initialise */
  currently_loaded_file[0]='\0';

  for (i=0;i<arguments.n_xml_files;i++){
    strncpy(filename,arguments.xml_files[i],BUFSIZE);
    printf("\nreading in XML report file %s\n",filename);
    if ((i==0)&&(arguments.output_filename[0]=='\0')){
      snprintf(autosave_filename,BUFSIZE,"%s.autosave",filename);
    } else if ((i==0)&&(arguments.output_filename[0]!='\0')){
      strncpy(autosave_filename,arguments.output_filename,BUFSIZE);
    }
    build_file_database(filename,&file_list);
  }
  if (arguments.inspect_transients==YES){
    printf("\nsearching for detected transient candidates\n");
    for (cycle_file=file_list;cycle_file!=NULL;
	 cycle_file=cycle_file->next){
      if (cycle_file->usable==YES){
	for (i=0,transient_count=0;i<cycle_file->n_peaks;i++){
	  if (cycle_file->transient_peaks[i]->transient_candidate==YES){
	    transient_count++;
	    if (transient_count==1){
	      printf("(%d) path: %s  file: %s  ",++nf,cycle_file->path,
		     cycle_file->filename);
	      if (cycle_file->usable==YES)
		printf("is usable\n");
	      else
		printf("is not usable\n");
	    }
	    /* 	    printf("  - observation time of transient %d = %ld\n",transient_count, */
	    /* 		   cycle_file->transient_peaks[i]->observation_time); */
	    strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",
		     gmtime(&(cycle_file->transient_peaks[i]->observation_time)));
	    printf("  - observation time of transient %d = %s\n",transient_count,
		   timestring);
	    printf("  - lag peaks at channels: ");
	    for (j=0;j<cycle_file->transient_peaks[i]->n_lags;j++){
	      printf("%d (%d,%d) ",
		     cycle_file->transient_peaks[i]->lag_peaks[j]->lag_channel,
		     cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_ant1,
		     cycle_file->transient_peaks[i]->lag_peaks[j]->baseline_ant2);
	    }
	    printf("\n");
	    /* call the menu routine to ask whether we should view this
	       transient peak */
	    option_made=transient_peak_menu(&cycle_file,i,YES);
	    if (option_made==-1){
	      /* reassign the peak as a non-transient candidate */
	      cycle_file->transient_peaks[i]->transient_candidate=NO;
	      xml_changed=YES;
	    }
	  }
	}
	if (transient_count>0){
	  printf("  + it has %d possible peaks, %d of which are transient candidates\n"
		 ,cycle_file->n_peaks,transient_count);
	}
	/* write out the XML file (autosave it), so we can always leave
	   and come back without having to look at the peaks we have already
	   looked at again */
	if (xml_changed==YES){
	  printf("autosaving XML file\n");
	  write_xml_file(autosave_filename,&file_list);
	  xml_changed=NO;
	}
      }
    }
  } 
  if (arguments.search_for_file==YES) {
    for (i=0;i<arguments.n_files_to_find;i++){
      cycle_file=file_info(arguments.files_to_find[i],&file_list);
      if (cycle_file!=NULL){
	printf("path: %s  file: %s  ",cycle_file->path,
	       cycle_file->filename);
	if (cycle_file->usable==YES)
	  printf("is usable\n");
	else
	  printf("is not usable\n");
	printf("  + it has %d possible peaks, %d of which are transient candidates\n"
	       ,cycle_file->n_peaks,transient_count);
      }
    }
  }
  if (arguments.make_file_list==YES){
    if ((list_of_files=fopen(autosave_filename,"w"))==NULL){
      printf("can't open %s for writing\n",autosave_filename);
    } else {
      for (cycle_file=file_list;cycle_file!=NULL;
	   cycle_file=cycle_file->next){
	/* go through each file in turn and output the name
	   of the file to the output file */
	if ((arguments.useful_files_only==NO)||
	    ((arguments.useful_files_only==YES)&&(cycle_file->usable==YES)))
	  fprintf(list_of_files,"%s/%s\n",cycle_file->path,
		  cycle_file->filename);
      }
      fclose(list_of_files);
    }
  }
  if (arguments.make_plots==YES){
    printf("\nmaking plots of potential transients\n");
    for (cycle_file=file_list;cycle_file!=NULL;
	 cycle_file=cycle_file->next){
      if (cycle_file->usable==YES){
	for (i=0,transient_count=0;i<cycle_file->n_peaks;i++){
	  if (cycle_file->transient_peaks[i]->transient_candidate==YES){
	    transient_count++;
	    /* we have found a transient candidate peak, so we make
	       a peak summary plot of it */
	    transient_peak_menu(&cycle_file,i,NO);
	  }
	}
      }
    }
  }

  cycle_file=file_list;
  while(cycle_file!=NULL){
    free_file=cycle_file;
    cycle_file=cycle_file->next;
    free_archive_file(&free_file);
  }

  exit(0);

}

int peak_summary_plot(correlation_data **all_correlations,int interactive,char *pgplot_device){
  float plot_minvp_x,plot_minvp_y,plot_maxvp_x,plot_maxvp_y;
  float minreal,maxreal,plot_row_height,plot_col_width;
  float straightline_x[2],straightline_y[2];
  int accept_plot=0,i,minchannel,maxchannel,xpanel,ypanel,ant1,ant2;
  int option_made;
  float *channels=NULL,baseline_length,old_height,x_string_length;
  float y_string_length;
  correlation_data *cycle_data=NULL;
  char baseline_string[BUFSIZE],timestring[BUFSIZE],title[BUFSIZE];
  char tmpfile[BUFSIZE];

  /* we assign the boundaries of the viewport */
  plot_minvp_x=0.05;
  plot_minvp_y=0.05;
  plot_maxvp_x=0.98;
  plot_maxvp_y=0.93;

  /* open the PGPLOT device */
  if (pgplot_device==NULL){
    /* by default we use the screen device */
    cpgopen("/xs");
  } else {
    cpgopen(pgplot_device);
  }
  cpgask(0);

  minchannel=0;
  maxchannel=(*all_correlations)->if_nfreq;
  MALLOC(channels,maxchannel*sizeof(int));
  for (i=minchannel;i<maxchannel;i++){
    channels[i]=i;
  }

  while(accept_plot==0){
    /* determine the axis limits, as every panel will
       have the same limits */
    cpgpage();
    minreal=maxreal=(*all_correlations)->pols_combined[0];
    for (cycle_data=*all_correlations;cycle_data!=NULL;
	 cycle_data=cycle_data->next){
      for (i=minchannel;i<maxchannel;i++){
	if (cycle_data->pols_combined[i]<minreal)
	  minreal=cycle_data->pols_combined[i];
	if (cycle_data->pols_combined[i]>maxreal)
	  maxreal=cycle_data->pols_combined[i];
      }
    }
    /* some plot parameters */
    plot_row_height=(plot_maxvp_y-plot_minvp_y)/PANEL_Y_ROWS;
    plot_col_width=(plot_maxvp_x-plot_minvp_x)/PANEL_X_COLS;
    /* now we go through the data and make the plots */
    for (cycle_data=*all_correlations;cycle_data!=NULL;
	 cycle_data=cycle_data->next){
      /* which panel are we going to plot this spectrum in? */
      which_panel(cycle_data->baseline,cycle_data->if_no,&xpanel,&ypanel);
      /* now set the viewport */
      cpgsvp(plot_minvp_x+xpanel*plot_col_width,
	     plot_minvp_x+(xpanel+1)*plot_col_width,
	     plot_minvp_y+ypanel*plot_row_height,
	     plot_minvp_y+(ypanel+1)*plot_row_height);
      /* and set the axis ranges */
      if (arguments.individual_limits==YES){
	minreal=maxreal=cycle_data->pols_combined[0];
	for (i=minchannel;i<maxchannel;i++){
	  if (cycle_data->pols_combined[i]<minreal)
	    minreal=cycle_data->pols_combined[i];
	  if (cycle_data->pols_combined[i]>maxreal)
	    maxreal=cycle_data->pols_combined[i];
	}
	cpgswin(minchannel,maxchannel,minreal,maxreal);
      }
      /* draw the axes */
      cpgsci(PGSCI_INVERSE);
      cpgsls(PGSLS_FULL);
      if ((xpanel==0)&(ypanel==0)){
	/* we want to label the amplitude and channel axes */
	cpgbox("BCNTS",0,0,"BCNTS",0,0);
      } else {
	cpgbox("BCTS",0,0,"BCTS",0,0);
      }
      /* draw the spectrum */
      cpgline((maxchannel-minchannel),channels,cycle_data->pols_combined);
      /* draw each individual polarisation in a different colour */
      cpgsci(PGSCI_RED);
      cpgline((maxchannel-minchannel),channels,cycle_data->individual_pols[0]);
      cpgsci(PGSCI_ORANGE);
      cpgline((maxchannel-minchannel),channels,cycle_data->individual_pols[1]);

      /* draw a vertical line along each peak channel */
      straightline_y[0]=minreal;
      straightline_y[1]=maxreal;
      if (interactive==YES){
	/* screen plot so light colour */
	cpgsci(PGSCI_YELLOW);
      } else {
	cpgsci(PGSCI_BLUE);
      }
      cpgsls(PGSLS_DASHED);
      for (i=0;i<cycle_data->spectrum_info.npeaks_found;i++){
	straightline_x[0]=straightline_x[1]=
	  cycle_data->spectrum_info.peak_channels[i];
	cpgline(2,straightline_x,straightline_y);
      }
      /* now put the baseline length in the panel */
      baseline_length=sqrtf(powf(cycle_data->u,2)+powf(cycle_data->v,2)+
			    powf(cycle_data->w,2));
      baseline_to_antennae(cycle_data->baseline,&ant1,&ant2);
      snprintf(baseline_string,BUFSIZE,"(%d,%d) %.2f m",ant1,ant2,
	       baseline_length);
      /* get the length of the baseline string so we know how to 
         position it */
      cpgqch(&old_height);
      cpgsch(0.6);
      cpglen(PGLEN_WORLD,baseline_string,&x_string_length,
	     &y_string_length);
      /* now we place the string at the top right of the viewport */
      cpgsci(PGSCI_INVERSE);
      cpgtext((float)maxchannel-1.2*x_string_length,
	      maxreal-(maxreal-minreal)/7.0,baseline_string);
      cpgsch(old_height);
    }
    /* put a title on the plot */
    strftime(timestring,BUFSIZE,"%0H:%0M:%0S %0d/%0m/%Y",
	     gmtime(&((*all_correlations)->observation_time)));
    strip_path(names_.file,tmpfile);
    snprintf(title,BUFSIZE,"%s %s @ %s",tmpfile,
	     (*all_correlations)->sourcename,timestring);
    cpgsvp(plot_minvp_x,plot_maxvp_x,plot_minvp_y,plot_maxvp_y);
    cpgswin(0,100,0,100);
    cpgsci(PGSCI_INVERSE);
    cpglab("","",title);

    if (interactive==YES){
      /* now we ask whether the user wants to change the plot */
      option_made=menu_question("(c)hange x-axis range, (a)ccept plot",
				2,1,'c','a');
      if (option_made==1){
	accept_plot=1;
      } else if (option_made==0){
	option_made=menu_question("lower x-axis value",-1,minchannel);
	minchannel=option_made;
	option_made=menu_question("upper x-axis value",-1,maxchannel);
	maxchannel=option_made;
      }
    } else {
      accept_plot=1;
    }
  }

  cpgclos();

  return(0);
}

int want_this_file(int *longest_baseline){
  /* this routine goes through the information in the RPFITS header and
     decides whether the file is the right type for our project.
     it will return 1 if the file is the right type, 0 otherwise */
  int i,retval=1,j,matched;
  char datobs[BUFSIZE],source_name[BUFSIZE],instrument[BUFSIZE],stokesname[BUFSIZE];
  time_t now_time,then_time,cutoff_time;
  double baseline_length,max_baseline_length=0;
  struct tm time_then,time_now;
  
  strncpy(instrument,names_.instrument,16);
  for (i=0;i<16;i++){
    if (instrument[i]==' '){
      instrument[i]='\0';
      break;
    }
  }
  instrument[15]='\0';
  strip_trailing_spaces(instrument);
  strncpy(datobs,names_.datobs,12);
  datobs[10]='\0';
  strip_trailing_spaces(datobs);
  for (i=0;i<if_.n_if;i++){
    for (j=0;j<if_.if_nstok[i];j++){
      strncpy(stokesname,names_.if_cstok+i*(if_.if_nstok[i])*2+j*2,2);
      stokesname[2]='\0';
    }
  }
  for (i=0;i<su_.n_su;i++){
    strncpy(source_name,names_.su_name+i*16,16);
    source_name[15]='\0';
    strip_trailing_spaces(source_name);
  }
  for (i=0;i<anten_.nant;i++){
    for (j=0;j<anten_.nant;j++){
      baseline_length=sqrt(pow((doubles_.x[i]-doubles_.x[j]),2)+
			   pow((doubles_.y[i]-doubles_.y[j]),2)+
			   pow((doubles_.z[i]-doubles_.z[j]),2));
      if (baseline_length>max_baseline_length){
	max_baseline_length=baseline_length;
	*longest_baseline=(256*anten_.ant_num[i])+anten_.ant_num[j];
      }
    }
  }
  
  /* OK, so we want files meeting the following criteria: */
  /* need to be using the ATCA */
  matched=0;
  for (i=0;i<CONSIDER_N_INSTRUMENTS;i++){
    if (strcmp(instrument,consider_instrument[i])==0)
      matched=1;
  }
  if (matched==0){
    retval=0;
  }
  
  /* need to have 2 IFs */
#ifdef MORE_IFS
  if (if_.n_if<MIN_N_IFS)
#else
  if (if_.n_if!=MIN_N_IFS)
#endif
  {
    retval=0;
  }
  
  /* these two IFs needed to be observed simultaneously */
#ifndef NOT_SIMULTANEOUS
  matched=1;
  for (i=1;i<if_.n_if;i++){
    if (if_.if_simul[0]!=if_.if_simul[i])
      matched=0;
  }
  if (matched==0){
    retval=0;
  }
#endif
  
  /* the observation date needs to be more than 18 months ago */
  time(&now_time); /* the current time */
  gmtime_r(&now_time,&time_now);
  /* we make a time_t that would be midnight UTC on the day AFTER the observation
     was taken */
  time_then.tm_sec=time_then.tm_min=time_then.tm_hour=0; /* midnight */
  sscanf(datobs,"%4d-%2d-%2d",&(time_then.tm_year),&(time_then.tm_mon),&(time_then.tm_mday));
  time_then.tm_year-=1900; /* years since 1900 */
  time_then.tm_mon-=1; /* months since january */
  time_then.tm_mday+=1; /* the day after the observation */
  then_time=timegm(&time_then);
  /* now we figure out 18 months before now */
  time_now.tm_mon-=18;
  cutoff_time=timegm(&time_now);
  /* now do the comparison */
  if (then_time>cutoff_time){
    retval=0;
  }
  
  for (i=0;i<if_.n_if;i++){
    /* need to have 32-channel data, ie. continuum observation */
#ifdef LESS_CHANNELS
    if (if_.if_nfreq[i]>(MAX_N_CHANNELS+1)) /* number of frequency channels is one less than that given */
#else
    if (if_.if_nfreq[i]!=(MAX_N_CHANNELS+1))
#endif
    {
      retval=0;
    }
    /* both IFs need to be at L-band */
    if (doubles_.if_freq[i]>MAX_IF_FREQUENCY){
      retval=0;
    }
    /* both IFs need to have a bandwidth greater than 64MHz */
    if (doubles_.if_bw[i]<MIN_BANDWIDTH){
      retval=0;
    }
  }

  return (retval);
  
}

int clear_correlations(correlation_data **all_correlations){
  /* this routine clears all the memory in the correlation list */
  correlation_data *free_correlation=NULL;

  while(*all_correlations!=NULL){
    free_correlation=*all_correlations;
    *all_correlations=(*all_correlations)->next;
    remove_correlation(&free_correlation);
  }

  return(0);
}

int remove_correlation(correlation_data **correlation_to_remove){
  /* this frees the memory used by a correlation element */
  int i;
  FREE((*correlation_to_remove)->vis);
  FREE((*correlation_to_remove)->wgt);
  for (i=0;i<(*correlation_to_remove)->if_nstok;i++){
    FREE((*correlation_to_remove)->stokes_name[i]);
    if ((*correlation_to_remove)->individual_pols!=NULL)
      FREE((*correlation_to_remove)->individual_pols[i]);
  }
  FREE((*correlation_to_remove)->stokes_name);
  FREE((*correlation_to_remove)->individual_pols);
  FREE((*correlation_to_remove)->pols_combined);
  FREE((*correlation_to_remove)->spectrum_info.peak_channels);
  FREE((*correlation_to_remove)->spectrum_info.peak_verified);
  FREE((*correlation_to_remove)->spectrum_info.peak_number);
  FREE(*correlation_to_remove);

  return(0);
}


int transient_peak_menu(archive_file **current_file,int current_peak,int interactive){
  int option_chosen,default_option,longest_baseline;
  int this_jstat,rpfits_result,baseline,flag,bin,if_no,sourceno;
  int keep_reading,wtf_result,j;
  float *vis=NULL,*wgt=NULL,ut,u,v,w;
  correlation_data *peak_correlations=NULL,*cycle_correlations=NULL;
  correlation_data *peak_transforms=NULL;
  char plot_filename[BUFSIZE];
  
  if (interactive==YES){
    option_chosen=menu_question("view this peak",2,1,'y','n');
    default_option=0;
  } else {
    option_chosen=0;
  }
  if (option_chosen==0){
    /* we want to view the transient peak */
    /* first we load the file */
    if (strlen((*current_file)->path)!=0){
      sprintf(names_.file,"%s/%s",(*current_file)->path,
	      (*current_file)->filename);
    } else {
      sprintf(names_.file,"%s",(*current_file)->filename);
    }
    if ((strcmp(names_.file,currently_loaded_file)!=0)&&
	(strcmp((*current_file)->filename,currently_loaded_file)!=0)){
      /* the currently loaded file is not the one we want, so 
         we have to open a new file and discard all the old data */
      clear_correlations(&current_file_correlations);
      current_file_correlations=NULL;
      this_jstat=JSTAT_OPENFILE;
      rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,
			      &flag,&bin,&if_no,&sourceno);
      if (this_jstat==JSTAT_UNSUCCESSFUL){
	/* can't open this file */
	printf("can't open file %s\n",names_.file);
	/* try to find this same file in the current directory */
	strcpy(names_.file,(*current_file)->filename);
	this_jstat=JSTAT_OPENFILE;
	rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,
				&flag,&bin,&if_no,&sourceno);
	if (this_jstat==JSTAT_UNSUCCESSFUL){
	  /* can't open this file */
	  printf("can't open file %s\n",names_.file);
	  return(-1); /* discard this peak */
	}
      }
      keep_reading=1;
      while(keep_reading){
	/* read the next header */
	this_jstat=JSTAT_READNEXTHEADER;
	rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,
				&flag,&bin,&if_no,&sourceno);
	if (this_jstat==JSTAT_SUCCESSFUL){
	  /* figure out if this data is useful */
	  wtf_result=want_this_file(&longest_baseline);
	  if (wtf_result==1){
	    /* yes it is useful */
	    read_data(&current_file_correlations,longest_baseline,&keep_reading);
	  }
	} else if (this_jstat==JSTAT_ENDOFFILE){
	  keep_reading=0;
	}
      }
      /* close the file */
      this_jstat=JSTAT_CLOSEFILE;
      rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,&flag,
			    &bin,&if_no,&sourceno);
      strcpy(currently_loaded_file,names_.file);
    }
    
    /* now we pick out the correlations that correspond to the 
       transient peak we want to look at */
    /* first get rid of any currently existing separate linkages */
    for (cycle_correlations=current_file_correlations;
	 cycle_correlations!=NULL;
	 cycle_correlations=cycle_correlations->next)
      cycle_correlations->separate_next=NULL;

    for (cycle_correlations=current_file_correlations;
	 cycle_correlations!=NULL;
	 cycle_correlations=cycle_correlations->next){
      if (cycle_correlations->observation_time==
	  (*current_file)->transient_peaks[current_peak]->observation_time){
	cycle_correlations->separate_next=peak_correlations;
	peak_correlations=cycle_correlations;
      }
    }

    /* compute the Fourier transforms for these data */
    change_to_lagspace(&peak_correlations,&peak_transforms);
    /* scalar add the polarisations */
    combine_polarisations(&peak_transforms,STOKES_I);

    /* we know where the peaks were from the log, so reinstate them */
    for (cycle_correlations=peak_transforms;cycle_correlations!=NULL;
	 cycle_correlations=cycle_correlations->next){
      for (j=0;j<(*current_file)->transient_peaks[current_peak]->n_lags;j++){
	if ((((*current_file)->transient_peaks[current_peak]->lag_peaks[j]->baseline_ant1*256)+
	     ((*current_file)->transient_peaks[current_peak]->lag_peaks[j]->baseline_ant2))==
	    cycle_correlations->baseline){
	  /* baselines match */
	  cycle_correlations->spectrum_info.npeaks_found=1;
	  MALLOC(cycle_correlations->spectrum_info.peak_channels,
		 cycle_correlations->spectrum_info.npeaks_found*sizeof(int));
	  cycle_correlations->spectrum_info.peak_channels[0]=
	    (*current_file)->transient_peaks[current_peak]->lag_peaks[j]->lag_channel;
	  break;
	}
      }
    }

    /* now plot the spectra */
    if (interactive==YES){
      peak_summary_plot(&peak_transforms,YES,NULL);
    } else {
      /* generate the name of the file to plot to */
      snprintf(plot_filename,BUFSIZE,"transient_summary_plot_%s_%s_peak%d.ps/cps",
	       (*current_file)->filename,
	       (*current_file)->transient_peaks[current_peak]->source_name,
	       (*current_file)->transient_peaks[current_peak]->peak_number);
      peak_summary_plot(&peak_transforms,NO,plot_filename);
    }
  }
  
  if (interactive==YES){
    option_chosen=menu_question("discard this peak",2,default_option,'y','n');
    if (option_chosen==0){
      /* tell the caller that this peak should be discarded */
      return(-1);
    } else if (option_chosen==1){
      /* tell the caller that this peak should be kept */
      return(0);
    }
  }

  return(0);
  
}

int read_data(correlation_data **all_correlations,int longest_baseline,int *keep_reading){
  /* this routine goes through an RPFITS file and reads in all the 
     data, putting it in a data structure for easy access later */
  int reading_data=-1,current_if=-1,this_jstat=-1,rpfits_result=-1;
  int baseline=-1,scannum=-1,flag=-1,bin=-1,if_no=-1,sourceno=-1;
  int correlations_added=0,add_result;
  float *vis=NULL,*wgt=NULL,ut=0,u=0,v=0,w=0;

  /* now we figure out how much data there is in this scan */
  reading_data=1;
  current_if=-1;
  scannum=0;
  while(reading_data==1){
    this_jstat=JSTAT_READDATA;
    rpfits_result=rpfitsin_(&this_jstat,vis,wgt,&baseline,&ut,&u,&v,&w,
			    &flag,&bin,&if_no,&sourceno);
    if (baseline==-1){
      /* calibration data, don't bother with it */
      /* but this does mean we need to change baselines */
      current_if++;
      if (current_if==if_.n_if){
	current_if=0;
	scannum++;
      }
      /* allocate the vis and wgt arrays to be the correct size */
      if (current_if!=-1){
	if (vis!=NULL){
	  FREE(vis);
	}
	if (wgt!=NULL){
	  FREE(wgt);
	}
	MALLOC(vis,((if_.if_nstok[current_if])*(if_.if_nfreq[current_if]))*2*sizeof(float));
	MALLOC(wgt,((if_.if_nstok[current_if])*(if_.if_nfreq[current_if]))*sizeof(float));
      }
    }
    if ((this_jstat==JSTAT_SUCCESSFUL)||
	(this_jstat==JSTAT_ENDOFSCAN)){
      if ((scannum>=SCAN_BOUNDARY_LIMIT)&&(baseline!=-1)){
	/* add this data to the structure */
	add_result=add_correlation(all_correlations,&vis,&wgt,baseline,ut,u,v,w,flag,bin,if_no,
				   sourceno,if_.if_nfreq[current_if],if_.if_nstok[current_if],
				   longest_baseline,scannum);
	if (add_result==0)
	  correlations_added++;
      }
    } else if (this_jstat==JSTAT_ENDOFFILE){
      reading_data=0;
      *keep_reading=0;
    } else
      reading_data=0;
  }
  FREE(vis);
  FREE(wgt);

  return(correlations_added);
}

int menu_question(char *question,int n_options,int default_option,...){
  /* this routine just asks a question and waits for a response from
     the user that matches the options supplied. the options should be
     given as single letters, eg. y,n,m,... the default option will
     be capitalised by this routine, and the position of the chosen
     option is returned */
  int option_chosen=-1,i;
  va_list ap;
  char *options=NULL,inputline[BUFSIZE];
  
  if ((n_options<1)&&(n_options!=-1))
    return(option_chosen);

  if (n_options>0){
    options=malloc(n_options*sizeof(char));
    
    printf("%s? ",question);
    printf("(");
    va_start(ap,default_option);
    for (i=0;i<n_options;i++){
      /* get the next option */
      options[i]=tolower(va_arg(ap,int));
      if (i!=0)
	printf("/");
      if (i!=default_option){
	printf("%c",options[i]);
      } else {
	printf("%c",toupper(options[i]));
      }
    }
    printf("): ");
    
    va_end(ap);
  } else if (n_options==-1){
    /* we accept a number from the user and return it */
    printf("%s? [%d] ",question,default_option);
  }

  /* wait for input */
  fgets(inputline,BUFSIZE,stdin);
  /* set for default option */
  option_chosen=default_option;
  if (n_options>0){
    /* only compare the first letter of the string */
    for (i=0;i<n_options;i++){
      if (tolower(inputline[0])==options[i]){
	option_chosen=i;
      }
    }
  } else if (n_options==-1){
    option_chosen=atoi(inputline);
  }

  /* return the chosen option */
  return(option_chosen);
}

int change_to_lagspace(correlation_data **all_correlations,
		       correlation_data **all_transforms){
  /* this routine takes the all_correlations list and make a new list which 
     will be the fourier transform of the all_correlations visibility */
  correlation_data *new_transform=NULL,*cycle_data=NULL;
  correlation_data *cycle_transforms=NULL,*second_cycle=NULL;
  fftw_complex *in=NULL,*out=NULL;
  fftw_plan p;
  int i,j;

  cycle_data=*all_correlations;
  if (cycle_data==NULL)
    return(0);

  /* prepare the FFT */
  if (cycle_data!=NULL){
    in=fftw_malloc(cycle_data->if_nfreq*sizeof(fftw_complex));
    out=fftw_malloc(cycle_data->if_nfreq*sizeof(fftw_complex));
    p=fftw_plan_dft_1d(cycle_data->if_nfreq,in,out,FFTW_BACKWARD,FFTW_ESTIMATE);
  }
  
  /* now go through the list */
  while(cycle_data!=NULL){
    /* prepare a new transform */
    MALLOC(new_transform,sizeof(correlation_data));

    /* copy the basic information over */
    copy_correlation_info(&cycle_data,&new_transform);
    /* copy the values to the in array for each polarisation */
    for (i=0;i<cycle_data->if_nstok;i++){
      for (j=0;j<cycle_data->if_nfreq;j++){
	in[j]=cycle_data->vis[j*2*cycle_data->if_nstok+i*2+0]+
	  cycle_data->vis[j*2*cycle_data->if_nstok+i*2+1]*I;
      }
      /* now do the FFT */
      fftw_execute(p);
      /* now copy the out array to the new structure */
      for (j=0;j<cycle_data->if_nfreq;j++){
	new_transform->vis[j*2*cycle_data->if_nstok+i*2+0]=crealf(out[j])/sqrt(cycle_data->if_nfreq);
	new_transform->vis[j*2*cycle_data->if_nstok+i*2+1]=cimagf(out[j])/sqrt(cycle_data->if_nfreq);
      }
    }
    /* add the new transform to the list */
    new_transform->next=NULL;
    new_transform->prev=NULL;
    new_transform->samebaseline_next=NULL;
    new_transform->samebaseline_prev=NULL;
    new_transform->separate_next=NULL;
    new_transform->raw_data=cycle_data;
    if (*all_transforms!=NULL){
      for (cycle_transforms=*all_transforms;cycle_transforms->next!=NULL;
	   cycle_transforms=cycle_transforms->next);
      cycle_transforms->next=new_transform;
      new_transform->prev=cycle_transforms;
    } else {
      *all_transforms=new_transform;
    }
    /* connect the baselines */
    for (second_cycle=new_transform->prev;second_cycle!=NULL;
	 second_cycle=second_cycle->prev){
      if ((second_cycle->baseline==new_transform->baseline)&&
	  (second_cycle->if_no==new_transform->if_no)){
	second_cycle->samebaseline_next=new_transform;
	new_transform->samebaseline_prev=second_cycle;
	break;
      }
    }
    cycle_data=cycle_data->separate_next;
  }

  /* deallocate the plan and free our memory */
  fftw_destroy_plan(p);
  fftw_free(in);
  in=NULL;
  fftw_free(out);
  out=NULL;

  return(0);
}

#define COPY_INFO_NUMBER(a,b,c) (*a)->c=(*b)->c
#define COPY_INFO_STRING(a,b,c) strcpy((*a)->c,(*b)->c)

void copy_correlation_info(correlation_data **original,
			   correlation_data **destination){
  /* this routine copies over the basic information (ie. everything but the vis array
     and the list pointers) from one correlation_data structure to another */
  int i,wgtsize;

  /* copy all the simple data */
  COPY_INFO_NUMBER(destination,original,baseline);
  COPY_INFO_NUMBER(destination,original,islongest_baseline);
  COPY_INFO_NUMBER(destination,original,ut);
  COPY_INFO_NUMBER(destination,original,observation_time);
  COPY_INFO_NUMBER(destination,original,u);
  COPY_INFO_NUMBER(destination,original,v);
  COPY_INFO_NUMBER(destination,original,w);
  COPY_INFO_NUMBER(destination,original,flag);
  COPY_INFO_NUMBER(destination,original,bin);
  COPY_INFO_NUMBER(destination,original,if_no);
  COPY_INFO_NUMBER(destination,original,sourceno);
  COPY_INFO_NUMBER(destination,original,if_nfreq);
  COPY_INFO_NUMBER(destination,original,if_nstok);
  COPY_INFO_NUMBER(destination,original,if_bandwidth);
  COPY_INFO_NUMBER(destination,original,scannum);
  COPY_INFO_STRING(destination,original,sourcename);

  /* copy the Stokes name data */
  MALLOC((*destination)->stokes_name,(*destination)->if_nstok*sizeof(char*));
  for (i=0;i<(*destination)->if_nstok;i++){
    MALLOC((*destination)->stokes_name[i],BUFSIZE*sizeof(char));
    COPY_INFO_STRING(destination,original,stokes_name[i]);
  }

  /* copy the weighting data */
  wgtsize=((*original)->if_nfreq) * ((*original)->if_nstok);
  (*destination)->wgt=NULL;
  (*destination)->vis=NULL;
  MALLOC((*destination)->wgt,wgtsize*sizeof(float));
  for (i=0;i<wgtsize;i++)
    COPY_INFO_NUMBER(destination,original,wgt[i]);
  /* allocate the vis array */
  MALLOC((*destination)->vis,2*wgtsize*sizeof(float));

  /* set some things to be NULL */
  (*destination)->individual_pols=NULL;

}

int combine_polarisations(correlation_data **all_correlations,int stokes_to_make){
  /* this routine goes through the list all_correlations and scalar adds the two
     polarisations together */
  correlation_data *cycle_data=NULL;
  int i,j;
  float average_amplitude,rms_noise;

  cycle_data=*all_correlations;
  if (cycle_data==NULL)
    return(0);

  while(cycle_data!=NULL){
    /* which Stokes parameter are we making */
    if (stokes_to_make==STOKES_I){
      /* Stokes I is just (XX + YY)/2 */
      /* allocate the necessary memory */
      CALLOC(cycle_data->pols_combined,cycle_data->if_nfreq,sizeof(float));
      MALLOC(cycle_data->individual_pols,cycle_data->if_nstok*sizeof(float*));
      for (i=0;i<cycle_data->if_nstok;i++){
	if ((strcmp(cycle_data->stokes_name[i],"XX")==0)||
	    (strcmp(cycle_data->stokes_name[i],"YY")==0)){
	  MALLOC(cycle_data->individual_pols[i],cycle_data->if_nfreq*sizeof(float));
	  for (j=0;j<cycle_data->if_nfreq;j++){
	    /* alright, we should note that if you combine the pols in this way after
	       doing the Fourier transforms, then actually you're not dealing with 
	       XX and YY anymore, so can't make Stokes I with (XX + YY)/2. Instead
	       what we are trying to do is reduce the amount of noise in the spectrum
	       by adding the two individual "polarisations" together, so we use the
	       square root of the squares */
/* 	    cycle_data->pols_combined[j]+=cycle_data->vis[j*2*cycle_data->if_nstok+i*2]/2; /\* we only take the real bit *\/ */
	    cycle_data->pols_combined[j]+=powf(cycle_data->vis[j*2*cycle_data->if_nstok+i*2],2); /* we only take the real bit */
	    cycle_data->individual_pols[i][j]=fabsf(cycle_data->vis[j*2*cycle_data->if_nstok+i*2]); /* we make it positive */
	  }
	}
      }
      strcpy(cycle_data->combined_pol_name,"I");
    }
    for (j=0;j<cycle_data->if_nfreq;j++)
      cycle_data->pols_combined[j]=sqrtf(cycle_data->pols_combined[j]);
    /* now we compute the RMS noise and average amplitude of the spectrum we just made */
    rms_noise=standard_deviation(&(cycle_data->pols_combined),
				 cycle_data->if_nfreq,&average_amplitude);
    cycle_data->spectrum_info.average_amplitude=average_amplitude;
    cycle_data->spectrum_info.rms_noise=rms_noise;
    cycle_data->spectrum_info.npeaks_found=0;
    cycle_data->spectrum_info.peak_channels=NULL;
    cycle_data->spectrum_info.peak_verified=NULL;
    cycle_data->spectrum_info.peak_number=NULL;

    cycle_data=cycle_data->next;
  }

  return(0);
}

void which_panel(int baseline,int if_no,int *xpanel,int *ypanel){
  /* this routine takes a baseline number and figures out which panel
     in the peak summary page it should be plotted in */
  int ant1,ant2,cx,cy,cif;
  /* the baselines are plotted like this: 
     12 13 14 15 16
     23 24 25 26 34  IF1
     35 36 45 46 56
     12 13 14 15 16
     23 24 25 26 34  IF2
     35 36 45 46 56
     xpanel=0, ypanel=0 is 35 IF2 (ie. bottom left)
     xpanel=4, ypanel=5 is 16 IF1 (ie. top right)
  */
  cx=0; cy=PANEL_Y_ROWS-1;
  for (cif=PANEL_MIN_IF;cif<=PANEL_MAX_IF;cif++){
    for (ant1=PANEL_MIN_ANT;ant1<PANEL_MAX_ANT;ant1++){
      for (ant2=ant1+1;ant2<=PANEL_MAX_ANT;ant2++){
	/* is this the baseline and IF we're after? */
	if ((baseline==(256*ant1+ant2))&&
	    (if_no==cif)){
	  /* yes it is */
	  *xpanel=cx;
	  *ypanel=cy;
	  return; /* no need to do any more */
	}
	cx++;
	if (cx==PANEL_X_COLS){
	  cx=0; cy--;
	}
      }
    }
  }
}

int add_correlation(correlation_data **correlation_list,float **vis,float **wgt,int baseline,
		    float ut,float u,float v,float w,int flag,int bin,int if_no,int sourceno,
		    int if_nfreq,int if_nstok,int longest_baseline,int scannum){
  /* this routine adds a correlation point to the correlation_data linked list at the
     end of the list */
  int i,nodata=0,wgt_size;
  correlation_data *cycle_data=NULL,*new_correlation=NULL;
  struct tm obs_time;
  char datobs[BUFSIZE];

  /* we don't add the correlation if the data is flagged as bad */
  if (flag!=0)
    return(-1);

  /* this loop gets us to the last entry in the list */
  if (*correlation_list==NULL){
    /* there is no data currently */
    nodata=1;
    cycle_data=NULL;
  } else 
    /* move cycle_data to the end of the list */
    for (cycle_data=*correlation_list;cycle_data->next!=NULL;cycle_data=cycle_data->next);
  
  /* allocate memory for our new entry */
  MALLOC(new_correlation,sizeof(correlation_data));
  wgt_size=if_nstok*if_nfreq;
  MALLOC(new_correlation->vis,2*wgt_size*sizeof(float));
  MALLOC(new_correlation->wgt,wgt_size*sizeof(float));

  /* initialise some pointers */
  new_correlation->spectrum_info.peak_channels=NULL;
  new_correlation->spectrum_info.peak_verified=NULL;
  new_correlation->spectrum_info.peak_number=NULL;
  new_correlation->individual_pols=NULL;

  /* now assign the values */
  /* the baseline number, and whether this is the longest baseline 
     in the array */
  new_correlation->baseline=baseline;
  if (baseline==longest_baseline)
    new_correlation->islongest_baseline=1;
  else
    new_correlation->islongest_baseline=0;

  /* the bandwidth of the observation */
  for (i=0;i<if_.n_if;i++){
    if (if_.if_num[i]==if_no){
      new_correlation->if_bandwidth=doubles_.if_bw[i];
      break;
    }
  }

  /* the time of the observation */
  new_correlation->ut=ut;
  obs_time.tm_sec=(int)ut; /* coz ut is given in seconds since midnight */
  obs_time.tm_min=obs_time.tm_hour=0;
  strncpy(datobs,names_.datobs,12);
  datobs[12]='\0';
  strip_trailing_spaces(datobs);
  sscanf(datobs,"%d-%2d-%2d",&(obs_time.tm_year),&(obs_time.tm_mon),&(obs_time.tm_mday));
  if (obs_time.tm_year>100){
    /* four digit year, need to subtract 1900 to fit into C's time model */
    obs_time.tm_year-=1900;
  } /* else the year is a two digit year, which is already the number of years since 1900 */
  obs_time.tm_mon-=1; /* months since January */
  new_correlation->observation_time=timegm(&obs_time); /* convert to simple time representation */
  
  /* the uvw plane location */
  new_correlation->u=u;
  new_correlation->v=v;
  new_correlation->w=w;

  /* is the data alright? */
  new_correlation->flag=flag;
  new_correlation->bin=bin;

  /* the IF number */
  new_correlation->if_no=if_no;

  /* details of the source being observed */
  new_correlation->sourceno=sourceno;
  for (i=0;i<su_.n_su;i++){
    if (sourceno==su_.su_num[i]){
      strncpy(new_correlation->sourcename,names_.su_name+i*16,16);
      new_correlation->sourcename[15]='\0';
      strip_trailing_spaces(new_correlation->sourcename);
      break;
    }
  }

  /* the names of the polarisations that were observed */
  new_correlation->stokes_name=NULL;
  MALLOC(new_correlation->stokes_name,if_nstok*sizeof(char*));
  for (i=0;i<if_nstok;i++){
    MALLOC(new_correlation->stokes_name[i],BUFSIZE*sizeof(char));
    strncpy(new_correlation->stokes_name[i],names_.if_cstok+(if_no-1)*if_nstok*2+i*2,2);
    new_correlation->stokes_name[i][2]='\0';
  }

  /* the number of frequency channels */
  new_correlation->if_nfreq=if_nfreq;

  /* the number of Stokes parameters observed */
  new_correlation->if_nstok=if_nstok;

  /* the number of the scan since the observation started */
  new_correlation->scannum=scannum;

  /* the data and the weight assigned to it */
  for (i=0;i<wgt_size;i++){
    new_correlation->wgt[i]=(*wgt)[i];
    new_correlation->vis[2*i]=(*vis)[2*i];
    new_correlation->vis[2*i+1]=(*vis)[2*i+1];
  }
  new_correlation->pols_combined=NULL;

  /* now attach it to the the end of the list */
  new_correlation->next=NULL;
  new_correlation->prev=cycle_data;
  new_correlation->separate_next=NULL;
  if (nodata==1)
    *correlation_list=new_correlation;
  else
    cycle_data->next=new_correlation;

  /* connect the baselines */
  for (cycle_data=new_correlation->prev;cycle_data!=NULL;
       cycle_data=cycle_data->prev){
    if ((cycle_data->baseline==new_correlation->baseline)&&
	(cycle_data->if_no==new_correlation->if_no)){
      cycle_data->samebaseline_next=new_correlation;
      new_correlation->samebaseline_prev=cycle_data;
      break;
    }
  }

  return(0);
}

float standard_deviation(float **values,int nelements,float *average){
  /* this routine calculates the standard deviation of the array passed to it */
  float value_stddev=0.0;
  int i;

  /* first get the mean value of the array */
  *average=meanvalue(values,nelements);

  /* now take the sum of each point minus the average */
  for (i=1;i<nelements;i++)
    /* DO NOT USE THIS ROUTINE AS A GENERAL PURPOSE STANDARD DEVIATION FINDER
       it ignores channel 0, since it will often be contaminated with a large peak
       whenever there is a continuum source at the phase centre of the observation */
    value_stddev+=powf(((*values)[i]-*average),2);
  /* divide by the number of elements */
  value_stddev/=nelements;
  /* and take the square root */
  value_stddev=sqrtf(value_stddev);

  return(value_stddev);
}

float meanvalue(float **values,int nelements){
  /* this routine calculates the average value of the array passed to it */
  float value_mean=0.0;
  int i;
  for (i=1;i<nelements;i++)
    /* DO NOT USE THIS ROUTINE AS A GENERAL PURPOSE AVERAGE FINDER
       it ignores channel 0, since it will often be contaminated with a large peak
       whenever there is a continuum source at the phase centre of the observation */
    value_mean+=(*values)[i];
  value_mean/=nelements;
  return(value_mean);
}

void baseline_to_antennae(int baseline_number,int *antenna1,int *antenna2){
  /* this routine takes a baseline number and gets the antenna numbers 
     from it, as baseline_number = (256*antenna1) + antenna2 */
  div_t result;
  result=div(baseline_number,256);
  *antenna1=result.quot;
  *antenna2=result.rem;
}

void strip_path(char *full_path,char *filename){
  /* this routine takes a full path name of a file and returns just the
     filename without the path information */
  int i;

  /* is there a slash in this filename? */
  if (strstr(full_path,"/")==NULL){
    /* no slash, so we just copy the whole string */
    strcpy(filename,full_path);
  } else {
    /* there is a slash */
    for (i=strlen(full_path)-1;i>=0;i--){
      if (full_path[i]=='/'){
	/* this is the last slash */
	strcpy(filename,full_path+i+1);
	break;
      }
    }
  }
}

/**
 * ConvertInput:
 * @in: string in a given encoding
 * @encoding: the encoding used
 *
 * Converts @in into UTF-8 for processing with libxml2 APIs
 *
 * Returns the converted UTF-8 string, or NULL in case of error.
 */
xmlChar *
ConvertInput(const char *in, const char *encoding)
{
    xmlChar *out;
    int ret;
    int size;
    int out_size;
    int temp;
    xmlCharEncodingHandlerPtr handler;

    if (in == 0)
        return 0;

    handler = xmlFindCharEncodingHandler(encoding);

    if (!handler) {
        printf("ConvertInput: no encoding handler found for '%s'\n",
               encoding ? encoding : "");
        return 0;
    }

    size = (int) strlen(in) + 1;
    out_size = size * 2 - 1;
    out = (unsigned char *) xmlMalloc((size_t) out_size);

    if (out != 0) {
        temp = size - 1;
        ret = handler->input(out, &out_size, (const xmlChar *) in, &temp);
        if ((ret < 0) || (temp - size + 1)) {
            if (ret < 0) {
                printf("ConvertInput: conversion wasn't successful.\n");
            } else {
                printf
                    ("ConvertInput: conversion wasn't successful. converted: %i octets.\n",
                     temp);
            }

            xmlFree(out);
            out = 0;
        } else {
            out = (unsigned char *) xmlRealloc(out, out_size + 1);
            out[out_size] = 0;  /*null terminating out */
        }
    } else {
        printf("ConvertInput: no mem\n");
    }

    return out;
}
